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Abstract 

In  this  paper,  we  develop,  analyze  and  test  a  local  discontinuous  Galerkin 
(LDG)  method  for  solving  the  Camassa-Holm  equation  which  contains  nonlinear 
high  order  derivatives.  The  LDG  method  has  the  flexibility  for  arbitrary  h  and 
p  adaptivity.  We  prove  the  L2  stability  for  general  solutions  and  give  a  detailed 
error  estimate  for  smooth  solutions,  and  provide  numerical  simulation  results  for 
different  types  of  solutions  of  the  nonlinear  Camassa-Holm  equation  to  illustrate 
the  accuracy  and  capability  of  the  LDG  method. 

AMS  subject  classification:  65M60,  35Q53 
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1  Introduction 

In  this  paper,  we  consider  numerical  approximations  to  the  Camassa-Holm  (CH)  equation 

[4,5] 


II /  'Ujxxt  +  2/S IUx  T  3uUx  ‘2'Ux'Ujxx  T  'U'^xxxi  (1*1) 

where  k  is  a  constant.  Such  nonlinearly  dispersive  partial  differential  equations  support 
peakon  solutions.  The  lack  of  smoothness  at  the  peak  of  the  peakon  introduces  high- 
frequency  dispersive  errors  into  the  calculation.  It  is  a  challenge  to  design  stable  and 
accurate  numerical  schemes  for  solving  this  equation. 

We  develop  a  class  of  local  discontinuous  Galerkin  (LDG)  methods  for  this  nonlinear 
CH  equation.  Our  proposed  scheme  is  high  order  accurate,  nonlinear  stable  and  flexible 
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for  arbitrary  h  and  p  adaptivity.  The  proof  of  the  L 2  stability  of  the  scheme  is  given  for 
general  solutions.  Error  estimates  are  given  for  smooth  solutions.  To  our  best  knowledge, 
this  is  the  first  provably  stable  finite  element  method  for  the  Camassa-Holm  equation. 

The  main  motivation  for  the  algorithm  discussed  in  this  paper  originates  from  the 
LDG  techniques  which  have  been  developed  for  convection  diffusion  equations  (contain¬ 
ing  second  derivatives)  [9]  and  nonlinear  wave  equations  with  high  order  derivatives 
[32,  33,  24,  27,  28,  29,  30].  In  these  papers,  stable  LDG  method  for  quite  general  nonlin¬ 
ear  wave  equations  including  multi- dimensional  and  system  cases  have  been  developed. 
The  proof  of  the  nonlinear  L 2  stability  of  these  methods  are  usually  given  and  success¬ 
ful  numerical  experiments  demonstrate  their  capability.  These  results  indicate  that  the 
LDG  method  is  a  good  tool  for  solving  nonlinear  equations  in  mathematical  physics. 

There  are  only  a  few  works  in  the  literature  for  error  estimates  of  the  LDG  method 
for  nonlinear  wave  equations  with  high  order  derivatives.  In  [31],  a  priori  error  estimates 
are  given  for  the  LDG  method  for  nonlinear  convection-diffusion  equations  and  KdV 
equations.  There  are  technical  difficulties  to  derive  the  L 2  a  priori  error  estimates  from 
the  cell  entropy  inequality  and  approximation  results,  because  of  the  possible  lack  of 
control  on  some  of  the  jump  terms  at  cell  boundaries,  which  appear  because  of  the 
discontinuous  nature  of  the  finite  element  space  for  the  DG  method.  The  remedy  in  [31] 
to  handle  such  jump  terms  is  via  a  special  projection,  which  eliminates  such  troublesome 
jump  terms  in  the  error  equation.  It  is  more  challenging  to  perform  L 2  a  priori  error 
estimates  for  nonlinear  PDEs  with  high  order  derivatives  than  for  first  order  hyperbolic 
PDEs  in  [34], 

The  CH  equation  (1.1)  was  derived  as  a  model  for  the  propagation  of  the  unidirec¬ 
tional  gravitational  waves  in  a  shallow  water  approximation,  with  u  representing  the 
free  surface  of  water  over  a  flat  bed  [4,  5,  18].  The  CH  equation  has  a  very  intriguing 
structure.  For  instance,  it  is  completely  integrable  and  models  wave  breaking  for  a  large 
class  of  initial  data.  This  equation  has  attracted  a  lot  of  attention  in  the  literature. 
Lenells  gave  a  detailed  discussion  of  its  conservation  law  in  [23]  and  also  classified  all  the 
weak  traveling  wave  solutions  in  [22],  Johnson  [19]  implemented  the  inverse-scattering 
transform  method  to  the  solution  of  the  CH  equation.  Li  and  Zhang  [25]  used  the  same 
method  to  solve  a  multiple-soliton  for  the  CH  equation.  Bressan  and  Constantin  [2,  3] 
developed  a  new  approach  for  its  analysis,  particularly  for  the  investigation  of  the  wave 
breaking. 

There  are  only  a  few  numerical  works  in  the  literature  to  solve  the  CH  equation. 
Holden  and  Raynaud  [17]  proved  the  convergence  of  a  finite  difference  method  for  the  CH 
equation  and  they  also  developed  a  convergent  numerical  scheme  based  on  multipeakon 
in  [16] .  Several  different  aspects  of  periodic  traveling  wave  solutions  of  the  CH  equation 
were  numerically  explored  in  [20]  and  the  error  analysis  of  a  spectral  projection  of  the 
periodic  CH  equation  was  given  in  [21],  A  finite- volume  method  was  developed  in  [1]  to 
simulate  the  dynamics  of  peakons. 

The  discontinuous  Galerkin  (DG)  method  we  discuss  in  this  paper  is  a  class  of  finite 
element  methods  using  completely  discontinuous  piecewise  polynomial  space  for  the  nu- 
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merical  solution  and  the  test  functions  in  the  spatial  variables.  The  DG  discretization 
results  in  an  extremely  local,  element  based  discretization,  which  is  beneficial  for  parallel 
computing  and  maintaining  high  order  accuracy  on  unstructured  meshes.  In  particular, 
DG  methods  are  well  suited  for  /ip-adaptation,  which  consists  of  local  mesh  refinement 
and/or  the  adjustment  of  the  polynomial  order  in  individual  elements.  They  also  have 
excellent  provable  nonlinear  stability.  The  LDG  method  for  the  Camassa-Holm  equation 
(1.1)  that  we  design  in  this  paper  shares  all  these  nice  properties.  More  general  informa¬ 
tion  about  DG  methods  can  be  found  in  [7,  10,  11,  12],  Recently,  Eskilsson  and  Sherwin 
[13,  14,  15]  also  presented  discontinuous  spectral  element  methods  for  simulating  ID 
linear  Boussinesq-type  equations,  dispersive  shallow  water  systems  and  2D  Boussinesq 
equations. 

The  paper  is  organized  as  follows.  In  Section  2,  we  present  and  analyze  our  LDG 
method  for  the  Camassa-Holm  equation  (1.1).  In  Section  2.2,  we  present  the  LDG 
method.  Details  related  to  the  implementation  of  the  method  are  described  in  Section 
2.3.  We  give  a  proof  of  the  L 2  stability  in  Section  3,  and  present  an  a  priori  error 
estimate  in  Section  4.  Section  5  contains  numerical  results  to  demonstrate  the  accuracy 
and  capability  of  the  methods.  Concluding  remarks  are  given  in  Section  6.  Some  of  the 
more  technical  proofs  of  several  lemmas  are  collected  in  the  Appendix. 

2  The  LDG  method  for  the  Camassa-Holm  equation 

2.1  Notation 

We  denote  the  mesh  by  Ij  =  [Xj_i,Xj+i],  for  j  =  1, . . . ,  N.  The  center  of  the  cell  is  x3  = 
+Xj+ 1)  and  the  mesh  size  is  denoted  by  h3  =  xJ+i  —  Xj_i,  with  h  =  maxi<j<jv  h3 
being  the  maximum  mesh  size.  We  assume  the  mesh  is  regular,  namely  the  ratio  between 
the  maximum  and  the  minimum  mesh  sizes  stays  bounded  during  mesh  refinements.  We 
define  the  piecewise-polynomial  space  14  as  the  space  of  polynomials  of  the  degree  up 
to  k  in  each  cell  Ij,  i.e. 

Vh  =  {v:  v  E  Pk{I3)  for  x  e  Ij,  j  =  l,...,N}. 

Note  that  functions  in  14  are  allowed  to  have  discontinuities  across  element  interfaces. 

The  solution  of  the  numerical  scheme  is  denoted  by  Uh,  which  belongs  to  the  finite 
element  space  14 .  We  denote  by  (■u/))+,i  and  (uh)~,  i  the  values  of  Uh  at  xj+i,  from 
the  right  cell  I3+i,  and  from  the  left  cell  Ij,  respectively.  We  use  the  usual  notations 
[uh]  =  —  %  and  Uh  =  k(uh  +  uh )  to  denote  the  jump  and  the  mean  of  the  function 

Uh  at  each  element  boundary  point,  respectively. 
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2.2  The  LDG  method 

In  this  section,  we  define  our  LDG  method  for  the  Camassa-Holm  equation  (1.1),  written 
in  the  following  form 


u  —  uxx  - 

= 

(2.1) 

Qt  +  f(u 

X  1 

)x  2^^)  XXX  2  ((Mx)  ) x 

(2.2) 

with  an  initial  condition 

u(x,  0)  =  Uo(x) 

(2.3) 

and  periodic  boundary  conditions 

u{ 

[x,  t )  =  u(x  +  L,t), 

(2.4) 

where  L  is  the  period  in  the  x  direction  and  f(u)  =  2 /vu+fw2.  Notice  that  the  assumption 
of  periodic  boundary  conditions  is  for  simplicity  only  and  is  not  essential:  the  method 
can  be  easily  designed  for  non-periodic  boundary  conditions. 

To  define  the  local  discontinuous  Galerkin  method,  we  further  rewrite  the  equation 
(2.1)  as  a  first  order  system: 

u  —  rx  =  q,  (2.5) 

r  —  ux  =  0. 


The  LDG  method  for  the  equations  (2.5),  where  q  is  assumed  known  and  we  would  want 
to  solve  for  u,  is  formulated  as  follows:  find  Uh,  rh  G  14  such  that,  for  all  test  functions 
p,  </>  G  14, 


uhpdx  +  /  rhpxdx  —  (rhp  )j+i  +  (?4p+)?_i  =  /  qhpdx, 


rh(j)dx  +  /  uh(j)xdx  -  ( uh(j)  )j+i  +  (uh(j)+)j_ 1  =  0. 


(2.6) 

(2.7) 


The  “hat”  terms  in  (2.6)-(2.7)  in  the  cell  boundary  terms  from  integration  by  parts  are 
the  so-called  “numerical  fluxes”,  which  are  single  valued  functions  defined  on  the  edges 
and  should  be  designed  based  on  different  guiding  principles  for  different  PDEs  to  ensure 
stability.  For  the  standard  elliptic  equation  (2.5),  we  can  take  the  simple  choices  such 
that 


rh  =  rh,  uh=u+,  (2.8) 

where  we  have  omitted  the  half-integer  indices  j+\  as  all  quantities  in  (2.8)  are  computed 
at  the  same  points  (i.e.  the  interfaces  between  the  cells).  We  remark  that  the  choice  for 
the  fluxes  (2.8)  is  not  unique.  We  can  for  example  also  choose  the  following  numerical 
flux 


4  =  r+,  uh  =  uh. 


(2.9) 
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For  the  equation  (2.2),  we  can  also  rewrite  it  into  a  first  order  system: 


Qt  +  /(«)*  -  Px  +  B(r)x  =  0, 

p-(b(r)u)x  =  0,  (2.10) 

r  —  ux  =  0, 


where  B(r)  =  |r2  and  b(r)  =  B'(r )  =  r.  Now  we  can  define  a  local  discontinuous 
Galerkin  method  to  equations  (2.10),  resulting  in  the  following  scheme:  find  qh,  Ph , 
Tfi  G  14  such  that,  for  all  test  functions  <p,  4)  rl  G  14, 


/  ( qh)t<fidx  -  /  (f(uh)  -ph  +  B(rh))<pxdx 
Jlj  Jl:i 

+  ((/ -  Ph  +  5(rft))4)i+i  -  ((/ -  +  S(rfc))^+)J-_i  =  0, 

[  Ph^dx  +  l  b(r h)uhipxdx  -  ( b(rh)uh 4  )i+i  +  (6(rh)ufcV’+)j-i  = 


rh(j)dx  +  /  uhr)xdx  -  (uhrj  )j+i  +  («h77+),-_i  =  0 


The  numerical  fluxes  in  equations  (2.11)-(2.13)  are  chosen  as 


Ph=Ph,  Sft  =  %)  =  B(rh),  b(rh)  = 


B(rt)  ~  ) 


o:  -  ri, 


,  uh  =  uj. 


(2.11) 

(2.12) 

(2.13) 


(2.14) 


Here  ()  is  a  monotone  flux  for  solving  conservation  laws,  i.e.  it  is  Lipschitz 

continuous  in  both  arguments,  consistent  (f(uh,Uh)  =  f(uh)),  non-decreasing  in  the 
first  argument  and  non- increasing  in  the  second  argument.  Examples  of  monotone  fluxes 
which  are  suitable  for  discontinuous  Galerkin  methods  can  be  found  in,  e.g.,  [8].  We  could 
for  example  use  the  simple  Lax- Friedrichs  flux 
^  _  l 

f(uh>uh)  =  2  (f(uh)  +  f(uh)-a(uh~uh))i  a  =  max  |  f(uh)l 


where  the  maximum  is  taken  over  a  relevant  range  of  Uh ■  This  Lax-Friedrichs  flux  is 
used  in  the  numerical  experiments  in  next  section.  The  definition  of  the  algorithm  is 
now  complete. 

We  remark  that  the  choice  for  the  fluxes  (2.14)  is  not  unique.  In  fact  the  crucial  part 
is  taking  ph  and  Uh  from  opposite  sides  and  B(rh )  and  Uh  from  opposite  sides. 


2.3  Algorithm  flowchart 

In  this  section,  we  give  details  related  to  the  implementation  of  the  method. 

•  First,  from  the  equations  (2.6)-(2.8),  we  obtain  q ^  in  the  following  matrix  form 

qh  =  A  uh,  (2.15) 

where  qh  and  Uh  are  the  vectors  containing  the  degrees  of  freedom  for  qh  and  Uh , 
respectively. 
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•  From  (2.11)-(2.14),  we  obtain  the  LDG  discretization  of  the  residual  —f{u)x  + 
\{u2)xxx  ~  |(( ux)2)x  in  the  following  vector  form 

(Qh)t  =  res  (uh).  (2.16) 

•  We  then  combine  (2.15)  and  (2.16)  to  obtain 

A  (uh)t  =  res  (uh).  (2.17) 

•  We  use  a  time  discretization  method  to  solve 

(uh)t  =  A_1res  (uh).  (2.18) 

This  step  involves  a  linear  solver  with  the  matrix  A.  We  perform  a  LU  decom¬ 
position  for  A  at  the  beginning  and  use  it  for  all  time  steps.  Any  standard  ODE 
solvers  can  be  used  here,  for  example  the  Runge-Kutta  methods. 

The  LDG  matrix  A  is  a  sparse  block  matrix,  hence  its  multiplication  with  vectors  and 
a  linear  solver  involving  it  as  the  coefficient  matrix  can  be  implemented  efficiently. 


3  L 2  stability  of  the  LDG  method 

In  this  section,  we  prove  the  L 2  stability  of  the  LDG  method  for  the  Camassa-Holm 
equation  defined  in  the  previous  section. 

Proposition  3.1.  (L2  stability)  The  solution  to  the  schemes  (2.6)-(2.8)  and  (2.11)-(2.14) 
satisfies  the  L2  stability 

L 

(uh  +  rl)dx<  0.  (3.1) 

Proof.  For  equation  (2.6),  we  first  take  the  time  derivative  and  get 

(rh)tpxdx  -  ((rh)tp~)j+ 1  +  ((rh)tp+)  ■  i  =  /  (< qh)tpdx .  (3.2) 

Since  (3.2),  (2.7)  and  (2.11)-(2.13)  holds  for  any  test  functions  in  14,  we  can  choose 


(uh)tpdx  + 


p  =  -uh,  f  =  (r  h)t,  <p  =  uh,  4  =  -rh,  rj  =  ph. 


With  these  choices  of  test  functions  and  summing  up  the  Eve  equations  in  (3.2),  (2.7) 
and  (2.11)— (2.13),  we  obtain 


{(uh)tUh  +  (■ rh)trh)dx  -  /  f{uh)(uh)xdx  +  (f(uh)u 


h  Jj+i 


(f(Uh)u 


h  )j-  = 


( PhUhfxdx  -  ( phuh  +  uhph  )j+!  +  (phUf  +  uhPh)j- 
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+ 


(B(rh)uh)xdx  +  (. B(rh)uh  +  b(rh)uhrh  )j+i  -  {B(rh)u+  +  b(rh)uhr^)j_ 
(( rh)tuh)xdx  -  {{rh)tui  +  uh(r^)t)j+i  +  {(rh)tu+  +  M/)(r+)t)i_ i  =  0. 


1 

2 


Taking  F(uh)  =  fUh  f(r)dr,  we  have 


/  (MtUh  +  (fh)trh)dx  +  T  +i  -  T  _i  +  0  _i  =  0, 

Jij  2  2  2 

where  the  numerical  entropy  fluxes  are  given  by 

^•+i  =  (~F(uh)  +  fuh  +Phuh  ~  (Ph^h  +  UhPh) 

~B(rh)uh  +  B(rh)Uh  +  b(rh)r~  +  (r^)tu^  -  (( rh)tu l  +  uh{r~)t\ 
and  the  extra  term  0  is  given  by 


(3.3) 


i+2 


0j_  1  =  (j-F(Wft)J  -  /[«/,]  -  [pftMfc]  +Ph[Mh]  +  Uh\ph] 

+  [B(rh)uh]  -  b(rh)uh[rh]  -  B(rh)[uh]  -  [(rh)tuh\  +  (rh)t[uh]  +  uh[(rh)t\ )  . 

With  the  definition  (2.8)  and  (2.14)  of  the  numerical  fluxes  and  after  some  algebraic 
manipulation,  we  easily  obtain 


-\Ph.Uh]  +  Ph[uh\  +  Uh[ph\  =  0, 

[B(rh)uh]  -  b(rh)uh[rh]  -  B(rh)[uh )  =  0, 

-[(rh)tuh]  +  ( rh)t[uh ]  +  uh[(rh)t]  =  0 

and  hence 

0,-1  =  ([CM]  -  /M),- *  >  o, 

where  the  last  inequality  follows  from  the  monotonicity  of  the  flux 
[p(uh)]  -  f[uh]  =  [  ( f(s )  -  f(ui,u£))ds  >  0. 

Juh 

Summing  up  the  cell  entropy  inequalities,  we  obtain  the  desired  L2  stability  (3.1).  □ 


4  Error  estimates  of  the  LDG  method 

In  this  section,  we  present  the  procedure  to  obtain  a  priori  error  estimates  of  the  LDG 
method  for  the  Camassa-Holm  equation. 
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4.1  Notations,  definitions  and  auxiliary  result 

In  this  section  we  introduce  notations  and  definitions  to  be  used  later  in  the  paper  and 
also  present  some  auxiliary  results.  We  first  review  an  important  quantity  measuring 
the  relationship  between  the  numerical  flux  and  physical  flux  introduced  in  [34,  31].  We 
then  define  some  projections  and  present  certain  interpolation  and  inverse  properties  for 
the  finite  element  spaces  that  will  be  used  in  the  error  analysis. 

4.1.1  Notations  for  different  constants 

We  will  adopt  the  following  convention  for  different  constants.  These  constants  may 
have  a  different  value  in  each  occurrence. 

We  will  denote  by  C  a  positive  constant  independent  of  h,  which  may  depend  on  the 
solution  of  the  problem  considered  in  this  paper.  Especially,  to  emphasize  the  nonlinear¬ 
ity  of  the  flux  f{u)  (or  other  nonlinear  fluxes),  we  will  denote  by  (7*  a  positive  constant 
depending  on  the  maximum  of  \  f"  |  or/and  |/w|.  we  remark  that  (7*  =  0  for  a  linear  flux 
f  =  cu  with  a  constant  c.  For  problems  considered  in  this  section,  the  exact  solution  is 
assumed  to  be  smooth  with  periodic  or  compactly  supported  boundary  condition.  Also, 
0  <  t  <  T  for  a  fixed  T .  Therefore,  the  exact  solution  is  always  bounded.  We  follow  the 
convention  [34]  to  redefine  the  nonlinear  functions  f(u),  B(u),  etc.  outside  their  ranges 
such  that  the  derivatives  of  these  nonlinear  functions  f'(u ),  f"(u ),  etc.  become  globally 
bounded  functions. 

4.1.2  A  quantity  related  to  the  numerical  flux 

For  notational  convenience  we  would  like  to  introduce  the  following  numerical  flux  related 
to  the  discontinuous  Galerkin  spatial  discretization.  f(co~,u+)  is  a  given  monotone 
numerical  flux  that  depends  on  the  two  values  of  the  function  cu  at  the  discontinuity 
point  Xj+ 1,  namely  cu^i  =  lu(x±+1).  The  numerical  flux  /(cu~,cu+)  satisfies  the  following 
conditions: 

(a)  it  is  locally  Lipschitz  continuous,  so  it  is  bounded  when  aA  are  in  a  bounded 
interval; 

(b)  it  is  consistent  with  the  flux  /(cu),  i.e.,  /(cu,cu)  =  /(cu); 

(c)  it  is  a  nondecreasing  function  of  its  first  argument,  and  a  nonincreasing  function 
of  its  second  argument. 

In  [34],  Zhang  and  Shu  introduced  an  important  quantity  to  measure  the  difference 
between  the  numerical  flux  and  the  physical  flux.  In  [31],  Xu  and  Shu  introduced  the  idea 
of  “uniform  dissipative  flux”.  For  completeness,  we  give  their  definition  in  the  following 
lemma. 


Lemma  4.1.  ^4l  For  any  piecewise  smooth  function  u  G  L2( 0, 1),  on  each  cell  boundary 
point  we  define 


a(f;uj)  =  a(f-u  ,ix+)  = 


M  VM  -  /M),  ifM  +  0; 

i\f(u)\,  if[w]  =  0, 


(4.1) 


where  f(ou)  =  /( u  ,ov+)  is  a  monotone  numerical  flux  consistent  with  the  given  flux  f. 
Then  a(f] u)  is  non-negative  and  bounded  for  any  (u~,u+)  G  M2.  Moreover  we  have 


\\fm 


<  a(/;w)  +  C„|[w]|, 

<  +Ct  |[oj]|2. 


(4.2) 

(4.3) 


Remark  4.1.  Examples  of  monotone  fluxes  which  are  suitable  for  the  discontinuous 
Galerkin  methods  can  be  found  in,  e.g.,  [8].  For  our  error  estimates,  we  rewrite  the 
numerical  flux  in  a  viscosity  form 

J(u~,uj+)  =  f(u~ )  +  f(u+)  -  A(cW,u;+)(w+  -  a;-)),  (4.4) 

and  assume  the  viscosity  coefficient  A(cu-,c<;+)  satisfies 


A(cu  ,  cu+)  >  Aq  >  0,  Aq  is  constant. 


(4.5) 


We  call  such  flux  as  a  “uniform  dissipative  flux”  [31].  The  well  known  Lax-Friedrichs 
flux  is  a  uniform  dissipative  flux  with  a  proper  choice  of  A.  This  property  is  necessary 
in  our  proof  because  of  a  lack  of  control  for  certain  jump  terms  at  cell  boundaries  due 
to  the  nonlinear  terms  and  the  high  order  derivative  terms. 

We  would  also  like  to  use  the  following  simplified  notation.  For  any  functions  c j  and 
0,  we  denote 

=  5^  <x(fiv)j+ i[0]2+i- 

1  <j<N 


4.1.3  Projection  and  interpolation  properties 

In  what  follows,  we  will  consider  the  standard  L2-projection  of  a  function  u  with  k  +  1 
continuous  derivatives  into  space  14,  denoted  by  V,  he.,  for  each  j, 

f  (Vu>(x)  —  uj{x))v(x)dx  —  0  Vu  G  Pk(Ij),  (4.6) 

Jij 

and  the  special  projection  V±  into  14,  which  satisfy,  for  each  j, 

(  (V+lo(x)  —  u(x))v(x)dx  =  0  Vu  G  Pk~1(Ij),  and  V+u(x+_  i )  =  u(x  _  i )  (4.7) 

4  3  * 

/  (V~u(x)  —  u(x))v(x)dx  =  0  Vv  G  Pfc_1(/j),  and  V~u{x~  fl)  —  oj(x,+  i). 

/ T  J+2  J  2 
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For  the  projections  mentioned  above,  it  is  easy  to  show  (c.f.  [6]) 

||u/||  +  hWu'Woo  +  h^\\uoe\\rh  <  Chk+\  (4.8) 

where  ue  =  Vu  —  u  or  ue  =  V±u  —  u.  The  positive  constant  C,  solely  depending  on  u, 
is  independent  of  h.  Th  denotes  the  set  of  boundary  points  of  all  elements  Ij. 

4.1.4  Inverse  Properties 

Finally  we  list  some  inverse  properties  of  the  finite  element  space  14,  that  will  be  used  in 
our  error  analysis.  For  any  Uh  G  14,  there  exists  a  positive  constant  C  independent  of 
Uh  and  h,  such  that 

(i)  1 1 <9*0^11  <  Ch~l\\ioh\\,  (ii)  \\uh\\Th  <  Ch~?\\uh\\,  (iii)  ||a;4oo  <  Ch~^\\ivh\\.  (4.9) 
For  more  details  of  these  inverse  properties,  we  refer  to  [6]. 

4.2  The  main  error  estimate  result 

We  state  the  main  error  estimates  of  the  semi-discrete  LDG  scheme  for  the  Camassa- 
Holm  equation. 

Theorem  4.2.  Let  u  be  the  exact  solution  of  the  problem  (2.1)-(2.2),  which  is  suffi¬ 
ciently  smooth  with  bounded  derivatives,  and  assume  f  G  C3 .  Let  Uh  be  the  numerical 
solution  of  the  semi-discrete  LDG  scheme  (2.6)-(2.8)  and  (2.11)-(2.1f)  and  denote  the 
corresponding  numerical  error  by  eu  =  u  —  Uh  and  er  —  r  —  rh  where  r  =  ux  is  defined 
by  (2.5).  For  regular  triangulations  of  I  =  (0,1),  if  the  finite  element  space  14  is  the 
piecewise  polynomials  of  degree  k  >  2,  then  for  small  enough  h  there  holds  the  following 
error  estimates 


\\u-uh\\2 +  \\r-rh\\2  <Ch2k,  (4.10) 

where  the  constant  C  depends  on  the  final  time  T,  k,  ||n||fc+i,  || r || fc+i  and  the  bounds 
on  the  derivatives  |/^|,  m  =  1,2,3.  Here  ||n||fc+i  and  || r || fc+i  ore  the  maximum  over 
0  <t  <  T  of  the  standard  Sobolev  k  +  1  norm  in  space. 

Remark  4.2.  Although  we  could  not  obtain  the  optimal  error  estimates  0(hk+1)  for 
u  due  to  some  extra  boundary  terms  arising  from  high  order  derivatives,  numerical 
examples  in  Section  5  verify  the  optimal  order  0{hk+l)  for  u.  For  the  solution  r^,  our 
numerical  results  indicate  that  k- th  order  accuracy  is  sharp. 

4.2.1  The  error  equation 

In  order  to  obtain  the  error  estimate  to  smooth  solutions  for  the  considered  semi-discrete 
LDG  scheme  (2.6)-(2.8)  and  (2. ll)-(2. 14) ,  we  need  to  first  obtain  the  error  equation. 
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Notice  that  the  scheme  (2.6)  and  (2.11)-(2.13)  is  also  satisfied  when  the  numerical 
solutions  are  replaced  by  the  exact  solutions.  We  then  obtain  the  error  equation 


J  (u-uh)tpdx-  J -  Qh)t(p  +  p)  +  (r  -  rh)(f>  +  (p-  ph )ip  +  (r  -  rh)rjjdx 
+  [  (r-rh)tpxdx-((rt-(rh)t)p-)j+i+((rt-(rh)t)p+)j_i 


+  /  (u  -  uh)(j)xdx  -  ((«  -  uh)(f>  )-,i  +  ((u-uh)(j)+),i 


+  (p~  Ph)Pxdx  -  ((p  -  ph)<p  ),,i  +  ((p  -  ph)<p+)j- 1 


+  (u-  uh)r]xdx  -  ((u  -  uh)p  )i+i  +  ((u  -  Uh)p+)  j-i 


/  (B{r)  -  B[rh))tpxdx  +  ((B(r)  -  B{rh))^)j+i  -  ((B(r)  -  5(r,))<^_± 

J  I 3 

+  /  (b(r)u  -  b(rh)uh)il>xdx  -  (( b(r)u  -  b(rh)uh)if>~)j+i  +  (( b(r)u  -  b(rh)uh)i/j+)j_ i 

Jr,  2  2 

-  f  (f(u)  -  f(uh))pXdx  +  ((/(«)  -  f)<p~)i+i  -  i(f(u)  -  7)^+),_i  =  o 


for  all  p,  4>,  ip,  i\),  rj  G  Vh. 
Define 


Bj(u  -  uh,  q-qh,p-  ph,  r  -  rh\  p,  (j),  <p,  ip,  rj) 

=  I  (u  -  uh)tpdx  -  J  ((q  -  qh)t(p  +  <p)  +  (r  ~  rh)<j)  +  (p  -  ph )ip  +  (r 

+  /  (r~rh)tpxdx-((rt-(rh)t)p-)j+i+((rt-(rh)t)p+)j_i 
■r  T:i 

+  [  (u  -  uh)(j)xdx  -  ((it  -  uh)<f>~)i+i  +  ((«  -  Uh)(p+)j-i 


(4.11) 


+  /  (p  ~  Ph)Pxdx  -  ({p  -  ph)(p  ),,i  +((p-ph)<p+)ii 


+  /  (u  -  uh)rjxdx  -  ((u  -  Uh)rj  )j+i  +  ((u  -  uh)r l+)j-±, 

J 


Hj(f]u,uh]<p)=  /  (f(u)  -  f(uh))tpTdx  —  ((f(u)~  f)tp  ).+  i  +  ((/(«)  -f)ip+)ji, 

J i, 

(4.12) 


B;r,u,rh,uh](p,ip) 


(4.13) 
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=  /  (B(r)  -  B(rh))<pxdx  -  ((B(r)  -  B(rh))<p~)j+ 1  +  ((5(r)  -  B(r,))/)3. 


(6(r)w  -  b(rh)uh)i/;xdx  +  ((6(r)«  -  b(rh)uh)ip  )j+±  -  ((b(r)u  -  b(rh)uh)^+):j_  i. 


Summing  over  j,  the  error  equation  becomes 

N 

Y  Bs(u  -  uht  q-qh,p~Ph,  r  -  )>, ;  p,  o.  ip,  e.  if) 


(4,14) 


1= 1 
AT 


=  J]  (^i(/;  <p)  +  R#> r>  uh ;  <a  0) 

3=1 

for  all  p,  0,  <p,  0,  Tj  G  14. 

Denoting 


S  =  "P+M  —  Uft,  se  =  V+u  —  u, 

£  =  Vq-qh,  ie  =  Vq-q,  (4.15) 

v  =  Pp  —  ph,  =  Vp  —  p, 
a  =  Vr  —  rh,  ae  =  Vr  —  r 

and  taking  the  test  functions 

P=  -s,  (f>  =  crt,  ip  =  s,  0  =  -a,  77  =  v, 

we  obtain  the  important  energy  equality 

N 

22  BAS  -  £  -  £*> v  ~ve,a  -  cre;  -s,  at,  s,  -a,  v)  (4.16) 

i=i 

N 

=  22  (^l(/i  w> Uh''  s)  +  Ri(&’  B;  r’ r/l’ s’  ~CT))  ’ 

l=i 


4.2.2  Proof  of  the  main  result 


In  this  subsection,  we  will  follow  the  idea  of  [31]  to  present  the  main  proof  of  Theorem 
4.2.  We  shall  prove  the  theorem  by  analyzing  each  term  of  the  energy  equation  (4.16). 

First,  we  consider  the  left  hand  side  of  the  energy  equation  (4.16).  The  estimates  for 
the  left  hand  side  of  the  energy  equation  are  given  in  the  lemma  below.  The  proof  of 
this  lemma  will  be  given  in  the  Appendix  A.l. 


Lemma  4.3.  The  following  equation  holds 

N 

22  ~  -  T,  v  -  ve,  a  -  <re;  -s,  at,  s ,  -a,  v) 

l=i 


N 


=  /  ( sst  +  ata)dx 


isdx  -22((ye  +  oet)[s\)j+\ 

i=i 


(4.17) 
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To  deal  with  the  nonlinearity  of  the  flux  f(u)  we  would  like  to  make  an  a  priori 
assumption  that,  for  small  enough  h,  there  holds 

|| u  —  uh  ||  <  h.  (4-18) 

We  will  justify  the  validity  of  this  a  priori  assumption  later.  For  the  linear  flux  f(u)  =  cu, 
this  a  priori  assumption  is  unnecessary. 

Corollary  4.4.  Suppose  that  the  interpolation  property  (4.8)  is  satisfied,  then  the  a 
priori  assumption  (4.18)  implies  that 

IKHoo  <  Chh  and  ||Qw  -  U/Jloo  <  ChS  (4.19) 


where  Q  =  V  or  Q  =  V±  is  the  projection  operator. 

Next,  we  consider  the  right  hand  side  of  the  energy  equation  (4.16).  We  can  rewrite 
it  into  the  following  form 


n  N  r 

J2Hj{f;u,uh]s)  =  Y^  /  (/ (u)  ~  f  (uh))sxdx 
j= i  j= i  Jlj 

N  N 

+  5^((/(w)  -  f{Uh))[s])j+  1  +  ^((/(Ufc)  -  f)[s])j+  1, 

3= 1  3= 1 

N 

E  B'>  r’  -<J) 

j=i 

AT  V 

+ 


(4.20) 


(4.21) 


i=i  o 


+  /  (&(*>  -  b{rh)uh)axdx  +  ^((6(  r)u-S(rh)u+)[a])i+i, 

i=i  Jlj  3= i 

where  we  take  into  account  the  periodic  boundary  condition  and  recall  the  average  Uh  is 
defined  by  Uh  =  \(ut  +  %)• 

The  estimates  for  the  equation  (4.20)  are  given  in  the  lemma  below. 

Lemma  4.5.  Suppose  that  the  interpolation  properties  (4-8)  are  satisfied,  then  we  have 
the  following  estimate  for  (4-20) 


(4.22) 


—  —  4°^’  +  (C  +  CUIMloo  4"  h  1 1  1 1  oo ) )  1 1  ^  1 1  ^  +  (C  +  C*h  1||eu||^0)h2fc+1. 

Remark  4.3.  For  the  proof  of  this  lemma,  we  refer  readers  to  Lemma  3.4  and  Lemma 
3.5  in  [31].  For  f(u)  =  2 nu  +  §w2  in  the  CH  equation  (1.1),  we  have  f"'{u)  =  0  and  we 
do  not  need  the  estimate  for  %  in  [31]. 
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The  estimates  for  the  equation  (4.21)  are  given  in  the  lemma  below.  The  proof  of 
this  lemma  will  be  given  in  the  Appendix  A. 2. 

Lemma  4.6.  Suppose  that  the  interpolation  properties  (4-8)  are  satisfied,  then  we  have 
the  following  estimate  for  (4-21) 

N 

i£*  j(b,  B ;  r,  u,  rh,  uh;  s ,  -a)  |  (4.23) 


3= 1 
N 


2fo+2 


<  ^(|&'(r><7c[<7]|  +  l6(r)(a6)  Nl )j+h  +  +  +  Ch 2 

1=1 

Now  we  are  ready  to  get  the  final  error  estimates  (4.10).  Combining  equations  (4.16), 
(4.17),  (4.22)  and  (4.23),  we  obtain 

r  1  i 


< 


/  (sts  +  <7t<7)da;  +  -«(/;  nh)[s]2 
Jo  4 

i-i  Jv  w 

/  stescb  +  ^((tTe  +  <7?)[s])i+i  +  J^(|&' 

”'0  7  =  1  7  =  1 


r)«cre[cr]|  +  |6(r)(cre)  [s]|)j+ 


j=i 

2 


+  (C  +  C*  ( 1 1  s  1 1  oo  +  h  1||eu||^0))||s||2  +  [C  +  (7*h  1||eu||^0)h"/l+1  +  C||(t||". 


Again  by  Young’s  inequality  and  the  interpolation  property  (4.8),  the  equation  becomes 

f1  1  ~ 

/  (sts  +  crt<j)dx  +  -a(f;uh) [s]2 


<(C  +  a(|NU  +  ^~i||e«||o0))llslr  +  (C1  +  C*h-l\\eu \Uh'2k+1  +  Ch2k  +  Cllall2. 

where  we  use  the  uniform  dissipation  property  of  the  numerical  flux  /.  Using  the  results 
(4.19)  implied  by  the  a  priori  assumption  (4.18)  and  the  positive  property  of  a(f)Uh), 
we  can  get  the  following  error  estimate 

r-l 


1  d 

2  dt  J0 


(. s 2  +  cr2)dx  <  C(||s||2  +  ||cr||2)  +  Ch: 


2k 


Thus  Theorem  4.2  follows  by  the  triangle  inequality  and  the  interpolating  property  (4.8). 

To  complete  the  proof,  let  us  verify  the  a  priori  assumption  (4.18).  For  k  >  2, 
we  can  consider  h  small  enough  so  that  Chk  <  J2h,  where  C  is  the  constant  in  (4.10) 
determined  by  the  final  time  T.  Then,  if  t*  =  sup{f  :  || u(t)  —  Uh(t)  ||  <  h},  we  would 
have  || u(t*)  —  Uh{t*)\\  =  h  by  continuity  if  t*  is  finite.  On  the  other  hand,  our  proof 
implies  that  (4.10)  holds  for  t  <  t*,  in  particular  || u(t*)  —  Uh(t*)\\  <  Chk  <  \h.  This  is 
a  contradiction  if  t*  <  T.  Hence  t*  >T  and  our  a  priori  assumption  (4.18)  is  justified. 

5  Numerical  results 


In  this  section  we  provide  numerical  examples  to  illustrate  the  accuracy  and  capability 
of  the  method.  Time  discretization  is  by  the  third  order  explicit  TVD  Runge-Kutta 
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method  in  [26].  This  is  not  the  most  efficient  method  for  the  time  discretization  to  our 
LDG  scheme.  However,  we  will  not  address  the  issue  of  time  discretization  efficiency 
in  this  paper.  We  have  verified  with  the  aid  of  successive  mesh  refinements,  that  in  all 
cases,  the  results  shown  are  numerically  convergent.  We  will  give  the  numerical  test 
results  for  the  CH  equation 

Ut  —  uxxt  +  3uux  =  2  uxuxx  +  uuxxx  (5.1) 

with  different  initial  conditions. 

Example  5.1.  Smooth  solution 

In  this  example,  we  test  the  scheme  with  smooth  traveling  waves.  Smooth  traveling 
waves  are  solution  of  the  form 

u(x,t)  —  4>{x  —  ct)  (5.2) 

where  0  is  solution  of  second-order  ordinary  differential  equation 

oi 

(ftxx  =  4*  T7  v7-  (5-3) 

VP  -  c) 

In  order  to  get  a  smooth  traveling  wave,  we  choose  the  constants  c  and  a  as  in  [16],  i.e. 
a  =  c  =  3.  The  initial  conditions  for  0  is 

m  =  1,  g(0)  =  0.  (5.4) 

It  gives  rise  to  a  smooth  traveling  wave  with  period  a  ~  6.46954603635.  We  use  fourth- 
order  explicit  Runge-Kutta  method  with  100000  points  to  approximate  the  solution  of 
the  equation  (5.3).  This  high  precision  solution  is  used  as  a  reference  solution  for  the 
smooth  traveling  wave.  The  L 2  and  L°°  errors  and  the  numerical  orders  of  accuracy  for 
u  at  time  t  =  0.5  with  uniform  meshes  are  contained  in  Table  5.1.  Periodic  boundary 
conditions  are  used.  We  can  see  that  the  method  with  Pk  elements  gives  a  uniform 
(AH-l)-th  order  of  accuracy  for  u  in  both  norms.  For  the  solution  r^,  the  accuracy  is  k-th 
order  in  both  norms  for  k  >  1  and  this  indicates  that  our  error  estimates  result  for  r  is 
sharp. 

Example  5.2.  Accuracy  test 

The  peakon  solutions  of  the  CH  equation  (5.1)  are  well  known  and  are  the  only 
traveling  waves  for  which  there  is  a  simple  explicit  formula.  The  peaked  traveling  wave 
solution  is 

u(x,  t)  —  ce~^x~a\  (5.5) 

where  c  is  the  wave  speed.  We  give  the  accuracy  test  with  c  =  0.25.  The  accuracy  is 
measured  in  smooth  parts  of  the  solution,  1/5  of  the  computational  domain  away  from 
the  peak.  The  L 2  and  L°°  errors  and  the  numerical  orders  of  accuracy  for  u  at  time  t  —  1 
with  uniform  meshes  in  [—25,  25]  are  contained  in  Table  5.2.  We  can  see  that  the  method 
with  Pk  elements  gives  a  uniform  (AH-l)-th  order  of  accuracy  for  u  in  both  norms. 
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Tabic  5.1:  Accuracy  test  for  the  CH  equation  (5.1)  with  the  exact  solution  (5.2).  Periodic 
boundary  condition.  Uniform  meshes  with  N  cells  at  time  t  =  0.5. 


u  - 

-  uh 

r  — 

Th 

N 

L2  error 

order 

L°o  error 

order 

L2  error 

order 

L°°  error 

order 

1.42E-01 

- 

3.08E-01 

- 

1.42E-01 

- 

3.08E-01 

- 

pO 

7.95E-02 

0.84 

1.77E-01 

0.80 

7.95E-02 

0.83 

1.77E-01 

0.57 

40 

4.23E-02 

0.91 

9.41E-01 

0.91 

4.23E-02 

0.94 

9.41E-02 

0.87 

80 

2.18E-02 

0.95 

4.83E-02 

0.96 

2.18E-02 

0.98 

4.83E-02 

0.97 

m 

1.16E-02 

- 

6.63E-02 

- 

1.16E-02 

- 

6.63E-02 

- 

p 1 

3.12E-03 

1.90 

1.86E-02 

1.84 

3.12E-03 

0.68 

1.86E-02 

0.24 

40 

8.05E-04 

1.95 

4.76E-03 

1.96 

8.05E-04 

0.85 

4.76E-03 

0.63 

80 

2.04E-04 

1.98 

1.19E-02 

2.00 

2.04E-04 

0.93 

1.19E-03 

0.87 

El 

1.41E-03 

- 

6.75E-03 

- 

1.41E-03 

- 

6.75E-03 

- 

p 2 

1.49E-04 

3.24 

9.06E-04 

2.90 

1.49E-04 

2.64 

9.06E-04 

2.64 

40 

1.70E-05 

3.13 

9.85E-05 

3.20 

1.70E-05 

2.06 

9.85E-05 

1.45 

50 

8.95E-06 

2.88 

4.96E-05 

3.07 

8.95E-06 

1.95 

4.96E-05 

1.77 

Table  5.2:  Accuracy  test  for  the  CH  equation  (5.1)  with  the  exact  solution  (5.5).  Periodic 
boundary  condition,  c  =  0.25.  Uniform  meshes  with  N  cells  at  time  t  —  1. 


N 

L 2  error 

order 

L°°  error 

order 

10 

8.71E-03 

- 

1.90E-02 

~ 

pO 

20 

2.18E-03 

2.00 

4.95E-03 

1.94 

40 

9.58E-04 

1.18 

2.36E-03 

1.07 

80 

4.08E-04 

1.23 

1.19E-03 

0.98 

10 

1.45E-02 

- 

3.17E-02 

- 

P 1 

20 

8.33E-04 

4.12 

2.06E-03 

3.94 

40 

1.14E-04 

2.87 

3.74E-04 

2.46 

80 

1.80E-05 

2.67 

8.82E-05 

2.08 

10 

1.60E-02 

- 

3.50E-02 

- 

p2 

20 

4.05E-04 

5.31 

1.23E-03 

4.83 

40 

2.81E-05 

3.85 

9.65E-05 

3.68 

80 

3.54E-06 

2.99 

1.29E-05 

2.90 

Example  5.3.  Peakon  solution 

In  this  example,  we  present  the  wave  propagation  of  the  periodized  version  of  the 
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solution  (5.5).  In  the  single  peak  case,  the  initial  condition  is 

c 


Uq(x)  =  { 


cosh (a/2) 
c 

[  cosh (a/2) 


cosh  (re  —  x0), 


x  —  x0\  <  a/2, 


(5.6) 


cosh(a  —  (re  —  re0)),  \x  —  rc0|  >  a/2, 


where  Xo  is  the  position  of  the  trough  and  a  is  the  period.  We  present  the  wave  propaga¬ 
tion  for  the  CH  equation  with  parameters  c  =  1,  a  =  30  and  xq  =  —5.  The  computational 
domain  is  [0,  a].  In  Figure  5.1,  the  peak  profile  at  t  —  0,  5,  10  and  the  space  time  graph 
of  the  solutions  up  to  t  =  10  are  shown.  The  lack  of  smoothness  at  the  peak  of  peakon 
introduces  high-frequency  dispersive  errors  into  the  calculation  and  will  cause  the  nu¬ 
merical  oscillation  near  the  peak.  In  our  computation  of  the  LDG  method,  we  use  the 
P5  element  with  N  =  320  cells  to  resolve  the  peak.  We  can  see  clearly  that  the  moving 
peak  profile  is  resolved  very  well. 


Example  5.4.  Two-peakon  interaction 

In  this  example  we  consider  the  two-Peakon  interaction  of  the  CH  equation  with  the 
initial  condition 


Uo(x')  =  <j>i(x)  +4>2(x), 


(5.7) 


where 


x  =  < 


cosh(a/2) 


cosh(a/2) 


cosh  (re  —  Xi), 


\x  —  Xi\  <  a/2, 


i  =  1,2.  (5.8) 


cosh(a.  —  (re  —  re*)),  |re  —  re/  >  a/2, 


The  parameters  are  c\  =  2,  C2  =  1,  x\  =  —5,  rc2  =  5,  a  =  30.  The  computational  domain 
is  [0,  a].  We  use  the  P 5  element  with  N  =  320  cells  in  our  computation  of  the  LDG 
method.  In  Figure  5.2,  the  two-peakon  interaction  at  t  —  0,  5,  12  and  18  are  shown.  We 
can  see  clearly  that  the  moving  peak  interaction  is  resolved  very  well. 


Example  5.5.  Three-peakon  interaction 

In  this  example  we  consider  the  three-Peakon  interaction  of  the  CH  equation  with 
the  initial  condition 


u0(x)  =  4>  i(x) +(t>2(x)  +  (f)3(x),  (5.9) 

where  (pi,  i  =  1,  2,  3  are  defined  as  in  (5.8).  The  parameters  are  Ci  =  2,  c2  =  1,  c3  =  0.8, 
aq  =  —5,  x2  =  —3,  rr3  =  —1,  a  =  30.  The  computational  domain  is  [0,  a].  We  use  the  P 5 
element  with  N  =  320  cells  in  our  computation  of  the  LDG  method.  In  Figure  5.3,  the 
three-peakon  interaction  at  t  —  0,  1,  2,  3,  4  and  6  are  shown.  We  can  see  clearly  that 
the  moving  peak  interaction  is  resolved  very  well. 
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t=5 


Figure  5.1:  The  peak  profile  of  the  CH  equation  (5.1)  with  the  initial  condition  (5.6). 
Periodic  boundary  condition  in  [0,30].  P5  elements  and  a  uniform  mesh  with  N  =  320 
cells. 

Example  5.6.  Solution  with  a  discontinuous  derivative 

In  this  example  we  consider  a  initial  data  function  uq  which  has  a  discontinuous 
derivative  as  in  [16].  The  initial  condition  is 

“oW  =  (sTPF  (5'10) 

The  computational  domain  is  [—30,30].  We  use  the  P 2  element  with  N  =  640  cells  in 
our  computation  of  the  LDG  method.  The  solutions  at  time  t  —  5,  10,  15  and  20  are 
shown  in  Figure  5.4.  Even  if  we  do  not  use  higher  order  polynomials  and  more  cells,  we 
still  obtain  a  good  resolution  of  the  solution  comparable  with  that  in  [16]. 

Example  5.7.  Break  up  of  the  plateau  traveling  wave 

In  this  example  we  consider  a  cut-off  peakon  in  [1],  i.e.  a  plateau  function  u(x,t )  = 
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Figure  5.2:  The  two-peakon  interaction  of  the  CH  equation  (5.1)  with  the  initial  condition 
(5.7).  Periodic  boundary  condition  in  [0,30].  P5  elements  and  a  uniform  mesh  with 
N  =  320  cells. 

4>(x  —  ct )  with 

(  cex+k,  x  <  —k, 

<f>(x)  —  <  c,  \x\  <  k,  (5-11) 

\  ce~x+k,  x  >  k. 

We  put  c  =  0.6  and  k  =  5.  The  computational  domain  is  [—40,40].  We  use  the  P2 
element  with  N  =  800  cells  in  our  computation  of  the  LDG  method.  The  break  up 
of  the  plateau  traveling  wave  is  shown  in  Figure  5.5  at  different  time.  The  solution  is 
resolved  very  well  comparing  with  the  result  in  [1], 


6  Conclusion 

We  have  developed  a  local  discontinuous  Galerkin  method  to  solve  the  Camassa-Holm 
equation.  L2  stability  is  proven  for  general  solutions,  and  an  a  priori  error  estimate 
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Figure  5.4:  The  solution  with  discontinuous  derivative  of  the  CH  equation  (5.1)  with 
the  initial  condition  (5.10).  Periodic  boundary  condition  in  [—30,30].  P 2  elements  and 
a  uniform  mesh  with  N  =  640  cells. 

is  obtained  for  smooth  solutions.  Numerical  examples  are  given  to  illustrate  the  accu¬ 
racy  and  capability  of  the  methods.  Although  not  addressed  in  this  paper,  the  LDG 
methods  are  flexible  for  general  geometry,  unstructured  meshes  and  h-p  adaptivity,  and 
have  excellent  parallel  efficiency.  The  LDG  method  has  a  good  potential  in  solving  the 
Camassa-Holm  equation  and  similar  nonlinear  equations  in  mathematical  physics. 

A  Appendix:  Proof  of  several  Lemmas 

A.l  Proof  of  Lemma  4.3 

Bj(s  -  se,  f  -  C,  v  ~  ve ,  a  -  ue;  -s,  at,  s,  -a,  v)  (A.l) 

=Bj (s,  f ,  v,  a;  -s,  at,  s,  -a,  v)  -  Bj(se,  £e,  ve,  ae;  -s,  au  s,  -a,  v), 
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.5:  The  bn 
mclition  (5. 
mesh  with 


By  the  same  argument  as  that  used  for  the  L 2  stability  in  Proposition  3.1,  the  first  term 
of  the  right  hand  side  in  (A.l)  becomes 


(■ sts  +  ata)dx  +  ^-+i  -  1 , 


where  T  =  v  s  —vs  —  sv  +  at  s  —  ats  —  'sat  . 

As  to  the  second  term  of  the  right  hand  side  in  (A.l),  we  have 


(A.2) 


Bj  (se,  C,  ve,  ae]  -s,  at,  s ,  -a,  v) 

=  /  (s^s  +  <7e(Jt)dx  +  /  (creu  —  vea)dx  +  /  (cr^ sx  +  se(at)x  +  vesx  +  sevx)dx  (A. 3) 

Jl;  Jlj  Jl; 

+  ((^e  +  <Ae)[s])i-i  + 


where  <h  =  —ves~  —  sev~  —  —  sea^ .  Because  V  is  a  local  L 2  projection,  and  P+, 

even  though  not  a  local  L2  projection,  does  have  the  property  that  s  —  V+s  is  locally 
orthogonal  to  all  polynomials  of  degree  up  to  k  —  1,  we  have 


aeatdx  +  /  (<rev  —  vea)dx  +  /  (afsx  +  se(at)x  +  vesx  +  sevx)dx  =  0. 


Noticing  the  special  interpolating  property  of  the  projection  V  ,  we  also  have 

=  0. 

The  equation  (A. 3)  then  becomes 

C,  ve,  ae;  -s,  at,  s,  -a,  v )  =  f  setsdx  +  ((rTe  +  (xte)[s]),,_i  +  $i+i  -  i, 


Combining  the  above  equation  with  (A.2),  summing  over  j  and  taking  into  account  the 
periodic  boundary  condition,  we  obtain  the  desired  equality  (4.17). 


A.2  Proof  of  Lemma  4.6 

The  proof  of  Lemma  4.6  is  similar  to  that  of  Lemma  3.5  in  [31].  The  main  difference  is 
that  we  take  fy  =  \{r^  +  ry  ) ,  and  ry  as  the  reference  values  of  the  functions  ry  and 
Uh  at  each  boundary  point.  For  the  nonlinear  terms  6(r)u  and  B(r),  we  use  the  following 
Taylor  expansions 

B(r )  —  B{rh )  —b(r)a  —  ^ b'(r)a 2  —  b(r)ae  +  b'  {r)aae  —  ^b'  {r)(ae)2 , 

5(r)  -  ^(?A )  -  T^VXO2  -  XX  «T  +  &,(r)u-(ae)-  -  ^5/(r)((ae)-)2, 

where  we  use  the  property  B'(r)  =  6(r)  and  6^m^(r)  =  0,  m>  2. 
b(r)u  -  b(rh)uh 
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—b'{r)ua  +  b(r)s  —  b'{r)sa  —  b'(r)ucre  —  b(r)se  +  b'(r)scre  +  b'(r)crse  —  b'(r)crese, 
b(r)u  —  b(rh)Wh 

—b'(r)ua  +  b(r)s+  —  6/(r)(xs+  —  b'(r)uae  —  6(r)(se)+ 

+  b'{r)s+ae  +  b'{r)a(se)+  —  b'  {r)ae  (se)+ . 

These  imply  the  following  representation 


R jib,  B]  r ,  u,  r h,  Uh ;  s,  —a )  —  <Si  +  S-2  +  £3  +  S4  +  S§ 


(A.4) 


where 


n  r  n 

S]  —  /  b'(r)uaaxdx  +  y^(^(r)tta-[cr])?-H 

.7=1  -'h  .7=1 


+£  b(r)(scr)xdx  +  ^^(6(r)(s  +  [a]+a  [s]))i+i, 


i=i  "  ^ 


52  =  “  9  (  51  /  &'Cr)(^2Ma;  +  ^(6,(r)(2s+a[a]  +  (a  )2[s]))j+i 


J  =  1 


AT  „  JV 

S3  =  -  ^2  /  b’{r)u(Jeaxdx  -  ^(b' {;r)uae[a})j 
7=1  •Rj  7=1 


/  6(r)(oresa;  +  sea;c)da;-^(6(r)((se)+[a] +  (ae)  [s])^ 


i=i  ■’ 


TV  .  TV 

b\r)ae(sa)xdx  +  ^(&'(r)(s+he[a]  +  (ae)"cr'[s]))i+i 
.i  1  Jlj  3= 1 

TV  „  AT 

+  5Z  /  &,(r)sWxcia;  +  ^(6/(r)a(se)  +  [a])j+i, 


j  =  l  " 


cS5  =  -  i  b' {r)((<je)2 sx  +  2aesecrx)dx 

N 

+  ^(6'(r)(2(s')+a'H  +  (K)-)2W))j+i 


will  be  estimated  separately  later. 

•  The  (Si  term. 

After  a  simple  integration  by  parts,  it  is  easy  to  obtain 


Si  =  --£  /  ( b'(r)u)xa2dx  -  ^  /  ( b(r))xsadx  <  C||o-||2  +  -||s||2.  (A. 5) 


j= 1  17  ^ 


i=i  17 
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•  The  S'2  term. 


After  a  simple  integration  by  parts,  it  is  easy  to  obtain 

S-2  =  -  ^2  /  ( b'(r))xsa2dx  =  0. 
1  j= i  ’W 


(A.6) 


The  S3  term. 

We  can  rewrite  S3  into  the  following  form 


iV  „  iV  „ 

/  (b'(r)u  —  b' (rj)uj)aeaxdx  —  /  br  (rj)ujaeaxdx 

i= i  ^Ij  j=i 

AT  „  AT 

-  /  (Kr)  +  + 


i=i 

N 


3=1  Jll 


N 


N 


~  ^(b'(r)uae[a})j+ i  -  ^(&(r)(<7e)  [s])i+i  -  ^(6(r)(se)  +  [n])i+i . 

.i  I  "  3  !  "  ./  I 

The  second  term,  the  fourth  term  and  the  last  term  in  the  above  equation  are  zero 
by  the  definition  of  the  special  projection.  Because  of  \b'(r)w  —  b'(rj)uj\  =  0(h) 
on  each  element  Ij,  then  by  the  inverse  property  (i)  in  (4.9),  together  with  the 
interpolation  property  (4.8),  the  first  term  in  the  above  equation  is  estimated  by 

Nr 

I  ^  /  ( b'(r)u  —  b'(rj)uj)aeaxdx\  <  C||<7e|| ||oj|  <  Cj|a||2  +  Ch2k+2. 

3  =  1  ’L‘ 

By  the  same  argument  we  can  also  get  the  estimate  for  the  third  term 

N  r  1 

I  ^  /  (6(r)  ~  &(ri))(°'e's^  +  seax)dx |  <  Cj|oj|2  +  -||s||2  +  Ch2k+2 . 

U  4  8 

Now  we  can  get  the  estimate  of  S3 

N  1 

<  ^(\b'(r)uv‘{*}\  +  |Kr)(«r')-[S]|)i+.  +  C||<r||2  +  -||S||2  +  Cfc2t+2.  (A.7) 

3=1 

Si  and  S5  terms. 


Because  S4  and  S5  are  high  order  terms  in  the  Taylor  expansion,  it  is  easy  to  show 
by  Young’s  inequality  and  the  inverse  properties  (i)  and  (ii)  in  (4.9)  that 


S4  <  Ch  1(||cre||oo  +  || "5e || 00) ||^"||  1 1 5 1 1  <  C||oj|2  +  -||s||2, 


(A.8) 


5,  <  Ch~1\\ae\ 


IMIIIXI  +  MINI)  <  c|M|2  +  i||s||2  +  ch2k+2.  (A. 9) 


Therefore,  summing  up  the  above  estimates  from  equations  (A. 5)  to  (A. 9),  we  complete 
the  proof  of  Lemma  4.6. 
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